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We consider an S — 1 Heisenberg chain with strong exchange (A) and single-ion uniaxial anisotropy (D) in 
a magnetic field (B) along the symmetry axis. The low energy spectrum is described by an effective 5=1/2 
XXZ model that acts on two different low energy sectors for a given window of fields. The vacuum of each 
sector exhibits Ising-like antiferromagnetic ordering that coexists with the finite spin stiffness obtained from the 
exact solution of the effective XXZ model. In this way, we demonstrate the existence of a spin supersolid phase. 
We also compute the full A — B quantum phase diagram by means of a quantum Monte Carlo simulation. 



PACS numbers: 75.10Jm, 75.40.Mg, 75.40.Cx 

One of the primary goals in the study of strongly corre- 
lated systems is the search for novel states of matter. These 
novel states usually emerge from competing interactions, each 
of which tends to stabilize distinct orderings. Different out- 
comes can result from such competition: a) none of the com- 
peting phases prevail, b) one interaction becomes dominant 
at low energies and stabilizes the corresponding ordering, c) 
more than one competing orderings coexist in a new phase. 
In the last case, the coexistence can be homogeneous or in- 
homogeneous. The supersolid phase is one of the simplest 
examples of homogeneous coexistence of two different phases 
(solid and superfluid) inbosonic systems. However, this novel 
state has been elusive both from an experimental [1J and a 
theoretical [2| viewpoint. While it is still unclear whether a 
supersolid (SS) state can be stabilized in the continuum, there 
are several numerical studies which show that a SS phase can 
be stabilized in the presence of a periodic potential or under- 
lying lattice EHUEJ. Surprisingly, this numerical effort has 
not been complemented whatsoever with rigorous analytical 
treatments that can demonstrate the existence of the SS phase, 
at least in certain limits. 

The SS state is easier to stabilize on a lattice because the 
lattice parameter of the "solid phase" or charge density wave 
cannot relax to any arbitrary value (it has to be an integer mul- 
tiple of the underlying lattice parameter). There are two nat- 
ural realizations of bosonic gases on a lattice: atomic gases 
in optical lattices or periodic substrates and quantum magnets 
(spin lattices). Although most of the previous numerical work 
has been focused on the former systems (3] 131 0, we have 
shown recently [6| that a class of spin Hamiltonians, which 
describe real magnets to a very good approximation, also con- 
tain a SS phase in their quantum phase diagram. In particular, 
the models for hard core bosons on frustrated lattices that have 
been considered in the context of atomic gases can also be 
generated as low energy effective theories of frustrated spin- 
dimer systems |6). The interest and the number of numerical 
works on spin SS phases is rapidly increasing OH] [9]]. 

The area of ID magnetism has also attracted a lot of atten- 
tion during the last few decades. This interest was mainly trig- 
gered by the successful application of field theory techniques 
for finding relevant and solvable low energy effective theo- 



ries. A remarkable accomplishment was achieved by Haldane 
ifTOll . who pointed out a qualitative difference between the 
low energy spectra of integer and half-odd-integer Heisen- 
berg chains. The gapped Haldane phase has peculiar proper- 
ties with measurable consequences, like the 5=1/2 end- 
chain states of open 5=1 Heisenberg chains ifTTI . 

In spite of the intense effort devoted to the study of ID 
Heisenberg like spin models, the search for new states and be- 
haviors is far from being complete. In this Letter, we demon- 
strate that the low energy spectrum of an S=l Heisenberg 
chain with uniaxial exchange and single-ion anisotropy con- 
sists of two sectors, each of which can be mapped into the 
exactly solvable (S=l/2) XXZ model for the limit of strong 
anisotropy. By exploiting this mapping, we also demonstrate 
that the ground state is a SS (Ising-like spin ordering coexist- 
ing with a finite spin stiffness) in a finite range of magnetic 
field. This is a surprising result considering that ID solid 
phases are usually unstable when the particle density (mag- 
netization for spin systems) deviates from the corresponding 
commensurate value. This occurs because each added (or re- 
moved) particle introduces a soliton in the ground state, some- 
thing that is particularly clear in the bosonization treatments 
lfT2ll . In our case, the solid (Ising) phase stabilized by the 
exchange anisotropy and the solitons remain as massive ex- 
citations for a finite range of densities due to the presence of 
a significant single ion anisotropy. This exceptional behavior 
leads to the stabilization of the supersolid state. In addition, 
we compute the full quantum phase diagram of H that cov- 
ers different regimes of anisotropy and was computed using a 
quantum Monte Carlo (QMC) algorithm. 

We start by considering an 5 = 1 Heisenberg model with 
uniaxial exchange and single-ion anisotropies on a linear 
chain C of L sites: 

H=Y, J±(S*Sf +1 +S«S? +1 )+J z S*S? +1 +DSf-BS*. 

i=l,L 

(1) 

Here J z > 0, D > and we define L + 1 = 1 to impose pe- 
riodic boundary conditions (PBC). Since H is invariant under 
spin-rotations along the z-axis, the z-component of the mag- 
netization, M = J2i is a good quantum number. We will 
assume that A = J z /J± > 1, Jj_ < D, J ± < J z - D and 
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J z > D. Therefore, we can treat J± as a perturbation. For 
J± = 0, H becomes diagonal in the basis of eigenstates of the 
set of operators {Sf }. The chain system consists of two inter- 
penetrating sublattices A and B. Since J z > D > 0, there are 
only two ground states for B = (i.e., M = 0): Sf — ±1 for 
i e A and Sf = Tl for i G 23: 



= ^LKIPn 00 o) 



where 1 00 0) is the reference state in which all the spins are 

in the eigenstate of Sf with eigenvalue zero: Sf jOO 0) = 

Vi G £. Any other state with M = has an energy 2J Z — D 
or higher. Let us consider now the case < M < L/2. In 
this case, the low energy subspace can again be divided into 
two different sectors or subspaces S A and S B . The subspace 
S A is generated by the following basis of states: 



2M/2 il 12 



(3) 



where the sites {ji..Jm} G S are all different: %\ < < 
iyi . In the same way, the subspace S B is generated by the 
basis of states: 



Q+ Q+ 

2M/2 ji ja' 



(4) 



where the sites {j"i...j'm} G .4 are all different: ji < j2--- < 
]m ■ The energy of any other state that has magnetization 
M and is orthogonal to S = S A © S B is higher by at 
least 2 J z — D. Therefore, we can use degenerate perturba- 
tion theory to solve the low energy spectrum for small J±: 
Jl < 2J Z - D and J ± < D. The subspaces S A and S B 
cannot be connected by any finite order process in the thermo- 
dynamic limit L — > oo. Consequently, we have two identical 
and disconnected low energy theories on each sector. From 
now on, we will consider one of them without loss of general- 
ity. The low energy subspace S A (Eq(3]) can be mapped into 
the Hilbert space for M hard core bosons on the L/2 sites of 
the B sublattice: 



b] b\ ....6J IO.4) 



(5) 



where i\ < i2---iN G B. The relation between the hard core 
bosons and the original spin operators is given by: 



s+ = V2b], sr = V2b { , s; = -i 



(6) 



where i G B and = b\b i . We emphasize that these rela- 
tions are only valid within the low energy subspace S A . In 
addition, Sf = 1 for i G A, which implies that the two-spin 
correlators (SfS?) and (S^~S~) have the following expres- 
sions up to quadratic corrections in the perturbative parameter 
J±/J z : 

(5|5| +r ) = e^(l - (nj) - (n J+r » + (n jnj+r ) il) 
(S+S- +r ) = 2{b]b ]+r ) (8) 



where b\ = and m = V« G A. After doing a canonical 
transformation and projecting out the high energy states: 

H = Pe- s He s P = H A + H B (9) 

we obtain the following expression for the low energy effec- 
tive model, Hj^, that acts on the sector S A : 

H A = LC + ^2 t(b]b l+2 + &l +2 6i) - mi + Vnim+2 (10) 



where L + 2 = 2 (PBC) and 
72 



c 



3 J z - 2D 



, t - 



J2 



2(J Z - D) 
72 

fi = -± -2t + AC + B + D - 2J Z , 

J 2 

V = -- ^ -2t + 2C*. 



(11) 



We note that \ V\ < —2t. Ha is the so-called t — V model (or 
S=l/2 XXZ Hamiltonian [14]) and is exactly solvable by the 
Bethe Ansatz method [ 13 1. The ground state is a Luttinger liq- 
uid for |V| < \2t\ and -2\t\ < fi < 2V+2\t\. The asymptotic 
behavior of the transverse and longitudinal two point correla- 
tors can be obtained with the bosonization method [ 12 1: 

h dr- 2K cos(27rpr) (12) 



3 j-\~ r 1 



p 2 



(b]b 



j+r/ 



= a 



2ir 2 r 2 
-g{K) 



< 2r -9(K) cog ar + c 3 r~ 1/2K cos nr (13) 
where g(K) = 2K + 1/(2K), a = [(2/3 - 1)tt , K is the 



Luttinger liquid (LL) parameter and p — 



Vi G B. A 



schematic contour map of K as a function of V/t and p can 
be found in Ref. lfT2l . The wave-length, A, of the oscilla- 
tions of the density-density correlator is the mean separation 
between bosons, i.e., A = l/p. Therefore, the longitudinal 
spin-spin correlator has two oscillatory components. The first 
component (first term of Eqj7j has a constant amplitude and 
the wave-length is equal to two lattice parameters. This con- 
tribution comes from the "solid" or Ising component of the 
SS phase. The second contribution (second term of Eq|7j) 
comes from the LL component and consequently decays with 
a power law according to Eq[T2] The density of the LL in 
the supersolid phase (p = p/2 since p the density on the B 
sublattice) can be extracted from the wave-length, 2A, of this 
oscillatory component. 

In order to test the accuracy of our effective model, we used 
the LANCZOS method to compute the exact ground state of 
the original Hamiltonian, H, in a L = 16 sites chain. Fig{T^ 
shows a comparison between the magnetization as function of 
field, M(B), obtained with the original and the effective mod- 
els for J z — 20 Jl and D = 10 Jl. The curves are practically 
indistinguishable. The full line is the M(B) curve obtained 
with Ha in the thermodynamic limit L — ► 00. FigJT]^ shows 
a similar comparison for the field dependence of the longitu- 
dinal spin structure factor: 



S zz (q) = -J2e- lqU ~ k) 



(SjS%.) 



(14) 
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Again, practically the same curves are obtained with H and 
Ha for a 16 sites chain. The finite value of S zz (q)/L for L — > 
oo indicates the existence of Ising ordering. Finally, FigfTJ; 
shows the same comparison for the field dependence of the 
spin stiffness p s or superfiuid density in the bosonic language. 
To define p s , we will use units of h = 1, lattice parameter 
a = 1 and effective mass of the bosons m = l/8t = 1. p s 
can be defined as the response to the gauge field generated by 
an infinitesimal flux <fi threading the ring (chain with PBC), 
Ps = §2^;, where e is the free energy per site (ground state 
energy density at T — 0). According to FigsJTJ? and[T];, the 
SS phase appears for B c \ < B < B c2 , where 



B cl = At + 2 J z - D - AC 



B r 



2V + 2J Z ~D- AC 



Jz 



(15) 



are determined by the conditions p c \ = 2t and p c2 — 2V—2t. 
Two different Ising phases appear on both sides: IS1 for B < 
B c i and IS2 for B c2 < B. In the dilute limit, B — > B^ (see 
Fig. [TJ;), the superfiuid density approaches the value of the 
total density p at T = 0. This is the known result for bosons 
in the continuum that is recovered in the dilute limit of the 
lattice system (pa — ► 0). 

The Ising order IS 2 corresponds to the state in which 
S? = 1 Vi € A and S? = Vi € B. If we keep in- 
creasing the field, ground states with higher values of M 
will be eventually stabilized. For L/2 < M < L, we can 
again identify a low energy subspace generated by the states: 
{\tPh..a m ) = ^S+S+....S+ M \00...0)} where i n denotes 
now any site of C and i\ < i^--- < iu- Again this subspace 
can be mapped into the one for hard-core bosons on the full 
chain C and the effective low-energy Hamiltonian is again a 
t — V model with t = J±, V = J z and p, = B — D. Since 
V ^> £, the IS2 phase (M — L/2) remains stable over a 
broad region of field. However, as shown in Fig(2] a second 
order transition to a spin liquid (SL) phase occurs at the criti- 
cal value B = B c ^ that is obtained from the exact solution of 
the Bethe-Ansatz equations for the t — V model ifTTl : 



B c3 = D + J ± A 



An Jj_ sinh 7 

7 



^ cosh \(2n + 1)tt 2 /2j} 
(16) 

where A = 2 cosh 7. The three transitions at B c \, B c a 
and £> c3 belong to the Dzhaparidze-Nersesyan-Pokrovsky- 
Talapov universality class |18|. Finally, the system becomes 
fully saturated at the critical field B C 4 = 2J Z + 2J± + D (this 
expression is also valid away from the strongly anisotropic 
limit A > 1). 

In order to extend the quantum phase diagram of H away 
from the strongly anisotropic limit, we have used the Stochas- 
tic Series expansion (SSE) [ 19] quantum Monte Carlo (QMC) 
method. The simulations were done on finite chains of length 
16 < L < 64. The SSE is a finite-temperature QMC based on 
the Taylor expansion of the partition function, e~P R . Ground 
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FIG. 1 : (Color online) Magnetization along the field direction M z 
(a), longitudinal structure factor S zz (n)/L (b), and spin stiffness p s 
(c) as a function of field for J z — 20 Jx and D — 10 J±. The 
dotted (dashed) line corresponds to the exact solution of the original 
(effective) model in a 16 sites chain. The full line corresponds to 
the exact solution of the effective model in the thermodynamic limit 
L — * 00. 



state estimates for the observables are obtained by choosing 
sufficiently large values of the inverse temperature [3. For 
the parameters explored in this study, (3 = 2L was found to 
be sufficient for the observables to have converged to their 
ground state values. The so-called Haldane state is the only 
new phase that appears in the full quantum phase diagram 
(see Fig j2j relative the strongly anisotropic limit. To char- 
acterize the different emergent phases, we computed M, p s 
and S zz (q). The spin stiffness, p s , is simply obtained by 
computing the winding number (W) fluctuations of the world 
lines:/9 s = (W 2 )/3.fM 

Both Ising phases, IS 1 and IS2, are marked by a finite value 
of S ZZ (Q) cx L and a vanishing value of p s in the limit 
L — > 00. The spin SS phase is characterized by a finite value 
of both S ZZ (Q)/L and p s in the same limit, while only p s re- 
mains finite (S ZZ (Q)/L — > 0) for the SL phase. Finally, the 
Haldane phase is characterized by a hidden ordering IF2T1 and 
both quantities, p s and S ZZ (Q) /L go to zero in the thermody- 
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namic limit. Since all these quantities are finite for finite size 
systems and estimates for L — > oo are obtained from finite- 
size scaling. 




FIG. 2: (Color online) Ground state phase diagram of the S — 1 
Heisenberg model in one dimension in the A — B parameter space 
with D — A/2. The spin-gapped Ising-like phases IS1 (m z = 
0) and IS2 (m z = 0.5) phases have long range diagonal order 
(S zz / L > for L — > oo) whereas the SL phase has a finite stiffness 
(p s > 0). The spin SS has simultaneous long range diagonal or- 
der and finite stiffness while both quantities are zero in the Haldane 
phase (HP). The symbols are results obtained from QMC simula- 
tions. The solid lines are the phase boundaries, B C i with 1 < i < 4 
(see Eqs|15|and|16} obtained from the effective low-energy models 
valid for A > 1. 

The results of the simulation are summarized in the ground 
state phase diagram for the A — B parameter space with 
D = A/2 (Fig(^). For A < 1, the ground state is a SL with fi- 
nite p s above a critical field B c q and up to the saturation field 
B C 4. Below B c o, the ground state is in the Haldane phase 
ifTOl . If we increase A, we find that the Haldane phase (HP) is 
separated from the IS1 phase by a line of critical points. This 
line appears because the excitations that becomes gapless at 
A = A c ~ 1.8 G2 and B = have M = (this is a second 
order transition between two M = ground states) while the 
\M\ > 1 excitations remain gapped. The Ising-like phase IS1 
evolves continuously into the state (or for A>1 
(see Eq|2]i. This corresponds to the empty band of the effec- 
tive model Ha. For A < 3.6, increasing the field B induces 
a transition from the IS 1 phase to the SL and from the SL to 
the fully polarized (FP) phase. At larger A(> 3.6), the Ising 
ordering persists along with the SL for a finite range of B giv- 
ing rise to a spin SS that is continuously connected with the 
SS phase that was analytically obtained in the limit A ^> 1. 
For A > 5.8, the SS phase ends up in the IS2 phase, as ob- 
tained for A ^> 1, while the SL phase appears in between 
for 3.6 < A < 5.8. At even higher fields, there is an IS2- 
SL transition at B = B c % and the system reaches saturation 
(M/L =1) for B = B c4 . Both phases and transitions are 
continuously connected with the ones obtained for A ^> 1. 
The phase boundaries, B C i with 1 < i < 4, obtained from 



the low energy effective models are shown in Fig |2] with solid 
lines. The close agreement with numerics for A ^> 1 confirms 
the validity of the low energy models in this limit. 

In summary, we have demonstrated the existence of a spin 
SS phase induced by field in a one-dimensional Heisenberg 
model with strong uniaxial anisotropy. This demonstration 
can be easily extended to higher dimensions We have also 
computed the quantum phase diagram as a function of the ex- 
change anisotropy A and field B. In particular, our ground 
state solution becomes asymptotically exact in the strongly 
anisotropic limit allowing for a full characterization of the SS 
spin phase. This is a remarkable result considering that SS 
phases found in other models were always obtained from nu- 
merical or approximated treatments. 
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